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Abstract — We performed accurate numerical calculations of angle-, time-, and frequency- dependent 
radiative transfer for the relativistic motion of matter in gamma-ray burst (GRB) models. Our technique 
for solving the transfer equation, which is based on the method of characteristics, can be applied to the 
motion of matter with a Lorentz factor up to 1000. The effect of synchrotron self-absorption is taken into 
account. We computed the spectra and light curves from electrons with a power-law energy distribution 
in an expanding relativistic shock and compare them with available analytic estimates. The behavior of the 
optical afterglows from GRB 990510 and GRB 000301c is discussed qualitatively. 
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INTRODUCTION 

The nature of the central sources of cosmic 
gamma-ray bursts (GRBs) has not yet been estab- 
lished. However, it is clear that GRBs with afterglows 
are at cosmological distances and release energy 
~10^^ erg on a time scale of the order of 100 s. 
The observed GRB peculiarities (nonthermal spectra, 
rapid temporal variability) require an ultrarelativis- 
tic motion of emitting plasma with characteristic 
Lorentz factors V ~ 100-300 (see Piran 2000; Blin- 
nikov 2000). 

In the standard GRB model (Rees and Meszaros 
1992), a photon— lepton fireball is produced (see 
Postnov 1998; Piran 2000). Initially, however, the 
GRB energy can also be electromagnetic (Usov 1994; 
Spruit 1999; Blandford 2002) and it probably propa- 
gates in a narrow cone (jet). The observed gamma- 
ray photons are generated by a nonthermal mecha- 
nism at the fronts of relativistic shocks (although the 
apparent nonthermal spectrum can also be explained 
in terms of the model of optically thick shells moving 
at relativistic velocities; see Blinnikov et al. 1999). 

Here, we develop a technique for solving the 
angle-, time-, and frequency-dependent transfer 
equation, which is based on the method of character- 
istics. It can be applied to the motion of matter with a 
Lorentz factor up to 1000. The main object of appli- 
cation of this technique must be the early generation 
phases of gamma-ray emission (during collisions 
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between internal shocks), for which the optical- 
depth effects can be noticeable. For now, however, 
we consider the radiation from the matter behind 
the front of an external shock and use an analytic 
solution (Blandford and McKee 1976) to describe the 
post-shock matter by taking into account the syn- 
chrotron self-absorption [cf. Downes et al. (2002), 
where the self-absorption was disregarded, but the 
hydrodynamics and spectrum of ultrarelativistic par- 
ticles were computed in a self-consistent way]. We 
computed the spectra and light curves from electrons 
with a power- law energy distribution in an expanding 
relativistic shock and compare them with available 
analytic estimates. 



FORMULATION OF THE PROBLEM 

One of the most popular models for GRB after- 
glows involves the propagation of a relativistic shell 
being decelerated by an external medium. The rela- 
tivistic shock heats up the captured matter as it enters 
the shell and causes the particles to be accelerated 
to ultrarelativistic energies. The X-ray and optical 
afterglows from GRBs in these models are associated 
with the nonthermal (synchrotron) radiation of rela- 
tivistic particles at the front of an external shock being 
decelerated in a circumstellar or interstellar medium 
(Meszaros and Rees 1997). Consider this problem 
in more detail by highlighting the most important 
points. 
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Fig. 1. The shape of the surface (the quasi-ellipsoid on the right) from which photons reach a remote observer (on the left) 
simultaneously. The explosion center is located at the vertex of the angle a. The farthest point of the visible surface lies at a 
small distance of ~(1 — (3) of the semimajor axis from the explosion center; the semiminor axis is ~l/7 of the semimajoraxis. 



The Propagation of Radiation from a Relativistic 
Shell 

Because of the high shock velocity, light from the 
ellipsoidal structure shown in Fig. 1 reaches the ob- 
server at a certain time. Let us determine the shape of 
the surface more accurately. 

Consider an emitting spherical shell of initial 
radius -R(to) = -Ro with an observer located at dis- 
tance D from its center. The shell begins to ex- 
pand as R = R{t). We assume that the time Iq 
at which the shell expansion begins corresponds 
to the time tgobs at which the observation begins; 
toobs = + {D - Rq) / c, where c is the speed of light. 

The radiation from points of a sphere with a radius 
depending on the cosine of the angle = cos a will 
reach the observer at some time tobs- For convenience, 
iobs is defined in such a way that it is equal to zero at 
the arrival time of the first signal on the shell motion. 
To determine the shape of the surface from which the 
radiation arrives, we take into account the fact that 
the time at which the propagation of photons from 
points of a sphere with a //-dependent radius R begins 
is the same and is specified only by tobs- 

In other words, 

"I — f^obs "T f^Oobs- 

C 

The surface shape can be determined from this 
equation by substituting in t = R^^{t). If we consider 
the simplest shell propagation equation 

R = Ro + pc{t -to), 0<P<1, 

where /? is the v/c ratio and 7 = (1 — /?^)^-^/^ is the 
Lorentz factor, then we obtain the equation of the 
surface 



R — R 

13c 



+ to + 



tobs "I" ^Oobs) 



R = 



As we see, in the approximation D R, this 
equation is the equation of an ellipse (Rees 1967). 
In Fig. I, the shape is more complex, because it 
corresponds to a variable velocity, as suggested by 
the solution of Blandford and McKee (1976). 

Each point on the sphere is characterized by the 
intensity Io{fi,r,UQ, cos 60) in the intrinsic frame of 
reference. In the observer's frame of reference, we 
write this intensity as /(/x, r, i^, cos 5) 

The intensity along the line of light propagation 
does not change in the absence of emission and ab- 
sorption sources and at the point of observation it 
will be the same as that at the point of radiation. 
Therefore, denoting /j,' = cos 9, we have for the flux 

1 

Fi, = 2tt J I{fi, r, u, cos S)fi'dfj,'. 

cos Smax 

Here, it is more convenient to pass from integration 
over to integration over a (Fig. 1 ): 

1 

F,^ = 2n J I{fj,, VyV, cos 5)fj,'{fj,)diJ,'{iJ,). 

/^min 

Denote R{n) /D = p{fi) and express cos 6 and d cos 6 
in terms ofp{n) and fi. The dependence R{fi) appears 
in the observer's frame of reference. It should be re- 
membered that p also depends on t and t, in turn, can 
be expressed in terms of iobs- However, to simplify 
our formulas, we will omit this dependence. For our 
subsequent calculations, we will need the following 
geometrical relations between the angles; 

cos 6*= — , \ dcos0 = p'^ 



cos 5 = 



ix--p{jj) 



COsJn 



cos (5-/3 
1 — /3 cos 5 ' 



{D > R). 



t^o 7(/^)(l - cos5/?(//))' 
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where /(//) = (1 +p^{n) - 2p{ii)iiY/'^. For the flux, 
we then have 

1 



obs j 



27r 



7 (i+p' 



//o (1) 



/^min 



X ^r(//), z^^— ^ , COS 5o(cos 5)j ^ — ^ d/Lt. 

When the flux is calculated, the condition imposed 
on the lower integration limit /Xmin is determined by 
the angle that corresponds to the maximum angular 
size of the shell from the point of observation: 

Our subsequent calculations are associated with 

a specific expression for the intensity /(r, i/q, cos(5o) 
on the shell surface and a specific shell propagation 
lawi?(t). 

The Transfer Equation 

For the intensity on the surface of a relativistic 
emitting shell to be calculated, we must solve the 
transfer equation in a comoving frame of reference. 
This is Eq. (2. 1 2 ) from Mihalas ( 1 980 ): 



(2) 



+ 7(1-m') 
-7'(/^ + /3) 



r c ot 

dl{fi,u) r/3(l-/x2) 



df3 



dr 



dfj, 



7 



7 



dp 



dp 



.2 



dl{fi, v) 



+ 37 



r c at 

dp 



Here, -q is the emission coefficient and x is the absorp- 
tion coefficient; the subscript was omitted, because 
all quantities refer to the comoving frame. 

Our numerical method of solution is described in 
the next section. 



Hydrodynamics 

The transfer equation (2) explicitly or implicitly 
(via ri and x) includes variables of the medium: its 
velocity, density, temperature, etc. For these vari- 
ables to be determined, we must solve the system 
of hydrodynamic equations. In general, the transfer 
and hydrodynamic equations constitute a combined 
system of equations. In our problem, however, we 



solve the transfer equation separately from the hydro- 
dynamic equations. At the same time, to determine 
the variables of the medium, we use a self-similar 
solution for a relativistic shock (with a Lorentz factor 
of the post-shock matter 7 » 1) in the spherically 
symmetric case for an ultrarelativistic gas (Blandford 
and McKee 1976). Let us give the formulas of this 
solution that we will need below. 

Taking the law of time variations in the shock- 
front Lorentz factor in the form oc t~"^ and choos- 
ing, for convenience, the self-similar variable C = [1 + 
2(m -I- l)r^](l — r/t), we derive the following expres- 
sions for the pressure, velocity, and density of the 
post-shock matter from the conditions at the shock 
front: 

P=lwir^f{0, 7' = k'5(C), (3) 



' = 2niT^h{C). 



n 

Here, wi is the enthalpy of the pre-shock matter, 
ni is its density, F is the shock-front Lorentz factor, 
and n' is the density of the post-shock matter in the 
observer's frame of reference. The matter density in 
the intrinsic frame of reference is related to the latter 
by n' = 7n. 

Substituting Eqs. (3) into the hydrodynamic 
equations yields a system of equations for /(C), ^(C), 
and /i(C) as a function of m. Consider the case m = 3, 
which corresponds to the conservation of total shell 
energy. The shock energy contained in the layer 
between the radii Ro{t) and Ri{t) is given by the 
expression 

Ri 

E{Ro,Ri,t)= J leirp'j'^r'^dr. 
Ro 

If we substitute in the solution for the functions /(C), 
5(C), and h{C) aim = 3: 

f = C'"'\ 9 = C\ h = c"\ 

then the total energy will be E = ^irwi x 
xt^r^/17, which gives the proportionality constant 
between F^ and t~^. 

Synchrotron Radiation 

An accurate calculation of the spectrum requires 
knowing not only the hydrodynamic quantities but 
also the electron energy spectrum and the magnetic- 
field strength. 

We assume that the electrons have a power-law 
distribution and that their total energy behind the 
shock front accounts for eg of the internal energy: 

iV(7) = Ko7~^,7 > 7mm,0 = 2 ' 
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where nie is the electron rest mass and Kq = {p — 

l)"07m7n,0- 

The magnetic field is parametrized by the quan- 
tity es, which is equal to the fraction of the inter- 
nal energy contained in the magnetic field: = 
STre^e. The magnetic field is randomly oriented and 
decreases with time due to the adiabatic shell ex- 
pansion. Other assumptions about the magnetic-field 
evolution and orientation weakly affect the resulting 
spectrum (Granot 1999). 

After the electrons have derived energy immedi- 
ately behind the shock front, they begin to lose it 
through adiabatic cooling determined by the solu- 
tion of Blandford and McKee (1976) and through 
synchrotron radiation. This process was described in 
more detail by Granot and Sari (2001). We present 
only the basic formulas for synchrotron radiation used 
in our calculations. 

The spectral power of a single electron averaged 
over the pitch angle is 



(-).. 



where 



R 



sy 



1 



CJTc5^(7e - 1 



Svr eB 



nipC 



and F{'u) is the standard function of synchrotron radi- 
ation (Rybicki and Lightman 1979). The synchrotron 
absorption coefficient is specified by the formula 



X 



^max 



NUMERICAL SOLUTION 
OF THE TRANSFER EQUATION 

The numerical solution of the problem is based on 
the simple and well-known method of characteristics 
(Mihalas 1980). We consider the relativistic transfer 
equation (2) in the spherically symmetric case in a 
comoving frame of reference. 

The main complexity of the equation is the pres- 
ence of four independent variables. The linearity of 
the equation allows its complexity to be decreased by 
constructing the characteristics for a given velocity 
field along which the differential operator is a total dif- 
ferential. If we choose the rays by describing them by 
a set of parameters and define s as some length along 
the ray, then we can determine the characteristics, the 
paths [t(s), r(s), /i(s), v{s)\, in such a way that 

dl drdi dadi dv dl dt dl 

— = h — 1 1 . 

ds ds dr ds dji ds dv ds dt 



We then derive the following system of equations that 
describe the characteristics from Eq. (2): 

dt 



7 



dfi 
ds 



ds 
dr 

ds 

2 



7(1 -M^) 



= 7(/^ + /3), 

l+^_7! 
r c 



di 



7'(m + /3) 



dp 



dv 
ds 



7 



7 



dr 

2 



dp 



+ 7V(/i + P) 



dp_ 
dr 



With the introduction of the characteristic rays, 
the transfer problem simplifies to 

dl{s) 



ds 



where 



X'(s) = X{s) + 37 



7y(s)-x'(s)/(s), 
Pil-li-") 



+ —{l + pfi)— + J fl{fi + p)— 



The characteristics are numerically computed by 

the fourth-order Runge— Kutta method with an adap- 
tive step. The step is adaptive, because at P close to 
unity, a small increment along the ray can lead to a 
significant angular jump comparable to the angular 
size of the emitting region. 

Note also that the variable quantities in the 
method must be of the order of unity, as follows from 
the constraint imposed on the cosine of the angle. 
Therefore, it is convenient to use a new system of 
units. We denote the physical and dimensionless 
(used in the program) quantities by the subscripts r 
and p, respectively. So, let 



Rr 



Tt 



nir 



Mnir, 



Denoting 

— Cp J I'f — Ip J 
Vr = JVp, Xr = XXp, 

where c is the speed of light and R, T, C, Y, J,and X 
are constants, we then derive the following expression 
for C, Y, J, and X in terms of R, T, and M after a 
simple work with dimensions: 

C = RT-^, Y = MT~^, 
J = MR-'^T-^, X = R-^. 
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The Analytic Solutions Used to Test the Numerical 

Method 

Below, we give some of the analytic solutions that 
we used to test the numerical method. 

Let us first consider the time- and frequency- 
independent transfer equation (2): 

7(M + /5)"^ + 7(l-M^) 



1 + 



dr 

-7'(/x + /3) 



+ 37 



+ 7V(M + /3)|^l/(r,/x) 



a/3 

dr 



r ■ - . 

V(.r,n) - x{r,tJ')I{r,fi). 



For constant emission and absorption coefficients 
and P = 0, the characteristics have the shape of 

straight lines r\f\^^]? = p (p is the parameter) and 
the equation has the analytic solution 

I(r,//) 



1 - exp { - X (/xr + Vi?^-r2(l-/x2)) } 



If we now assume /? to be constant, then the shape of 
the characteristics changes: p(l + = r\J\ — ji^, 
and the solution becomes slightly more complex, 

/(s) = /(r,M) 

7^r l^{p)R 



= -|l- 

where 



exp 



-X 



V/3V + (^2 + /3p2)(i22_p2 



]}. 



B? + /3V 

In the absence of absorption and for a constant 
emission coefficient on the right-hand side of the 
transfer equation, the analytic solution for the central 
characteristic (/x = ±1) is 



at/3: 



7(1-/3)^ 
const and 



/(O) = r]R 
I{R) = 27?i?(^ 



,ln(l + B) 



B 



1 



-S\3/2 



1 , l + i? ; 



where S = /3(i?), at /3 ~ r. 

To test the temporal component, we used the 
problem of intensity change near the surface of 
an emitting static transparent spherical shell with 
time when the emission coefficient abruptly changes 
from ?7o to -qi. If we use the following notation: p = 




Fig. 2. Instantaneous afterglow spectra at various times 
t = lO'^ s, where N is the number near the curve; Fm and 
Va are the analytically estimated flux and self-absorption 
frequency, respectively. 



R/D is the ratio of the shell radius to the distance 
from the shell center to the point of observation, 
t+ is the time measured from the beginning of the 
intensity change at the point of observation, c is the 
speed of light, and r = t+c/D, d = 1 — p + r, then 
the solution for this problem is 

i^i(r) = Gi(r?i,(l-p)2)-Gi(r?i,d2) 

+ Gi(%,d2)-Gi(%,l-p'), 



Gi{r],x] 



7rr,Ri2 {l-py {1 - p^f 
3 x3/2 xV2 



4p 
+ 2(1- 



p'^)x 



1/2 



r.3/2 



Here, r changes from to r* = ^/l — p"^ - 1 + i.e., 
within the interval during which the changes at the 
point of observation occur, and -Fi(r) is the flux. In 
the static case, this flux is 



p 

-1% 



-p){l- ^ip)p'^ 
+ p^ — 2p/x)2 



I{li*)dix, 



where 



fi-p 



(1 +p2 _ 2p//)V2 

is the cosine of the angle to the normal to the sphere 
surface on which the intensity depends. 

Comparison of the solutions considered above 
with the calculations based on our numerical method 
shows that the numerical method is applicable to the 
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lgt,s 

Fig. 3. Afterglow light curves for a set of frequencies 
If = lo" Hz, where N is the number near the curve. 

motion of matter up to Lorentz factors 7 ~ 1000 with 
an error less than 1%. 

RESULTS OF THE NUMERICAL SOLUTION 

We chose the following parameters for our nu- 
merical calculations of the afterglow spectra. These 
include the parameters that describe the hydrody- 
namics: the energy £"0 released through the process 
that leads to a GRB and the ambient density ni, the 
parameters that describe the radiation: the fraction of 
the internal energy contained in the magnetic field eb, 
the fraction of the internal energy transferred to elec- 
trons Ce, and the power-law index in the electron 
energy distribution p; and one more parameter: the 
photometric distance to the GRB D. 

Our main calculation, whose results are presented 
here, is based on the following published parame- 
ters: Eq = 10^^ erg, ni = 1 cm~^, eg = 0.5, eb = 0.1, 
p = 2.5, and D = lO^^^ cm. 

The large amount of released energy Eq is related 
to the spherical symmetry of the problem, while the 
observed GRBs can represent a jet with a solid an- 
gle J7. The total energy will then be lower by a factor 

The computed spectra and light curves are shown 
in Figs. 2 and 3. Here and below, the time is measured 
in the observer's frame of reference. Let us compare 
our results with available theoretical estimates (Hur- 
ley et al. 2002). In these estimates, the synchrotron 
spectrum is described by the maximum flux Fmax and 
three characteristic frequencies (fmin, t'c, z^a), where 
fmin is the synchrotron frequency of the electron with 
minimum energy whose Lorentz factor is 7niin,o, i^c is 
the cooling frequency, and i^a is the self-absorption 
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Fig. 4. Comparison of the frequencies that correspond 

to the maximum flux of the instantaneous spectra in the 
numerical calculations and analytic estimates. 

frequency. For these four parameters, the theoretical 
estimates are given by the formulas 

Ua = 2x KfHzEll^n/^e'^e^J^ = 4 x lO'^Hz, 

j/cool = 9 X lOi^HzEg-'/'nr'e-'/'t-y' 

= 2.66 X lO^^Hzi^^/^ 

l^mm — X iU nZ-C/52 ^e^B ''day 

= 3.80 X lO^^HzC^''^ 
Fmax = 20mJyE52ny'^e^^'^d^^ = 6.32 x lO-^^mJy. 

Between these frequencies, the spectrum is a 
power law with indices (2, 1.3, — 1/2, — p/2) for t < 
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Fig. 5. Comparison of the times that correspond to the 
maximum flux of the light curves in the numerical calcu- 
lations and analytic estimates. 
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r rg. b. 1 he intensity on the radiation suriace tor a set 
frequencies v = 10^ Hz, where N is the number near the 
curve. 

to = 4.2x10^ s and (2, 1.3, -(p - l)/2, -p/2) for 
t > to. 

Let us also compare the computed light curves 
with theoretical estimates (Sari et al. 1998), in which 

the characteristic times (tmirnicool) were calculated 
from the flux Fmax and the characteristic frequencies 

tcool = 7.3 X 10~^E^^^n^^e]^^iyfr^ day = GSu'^^ s, 
tniin = OmEli'et/hfiy-,^^' day 
= 2.37 X lO^iy^^^^' s. 

By introducing the frequency vq = fcooi(*o) = 
imin(io) = 1.14 X 10^^ Hz, we separate two cases: 
imin < icool for 1^ > vo and t^m > icooi for v < i^o- 

The results of our comparison are presented in 
Figs. 4 and 5. These figures show the frequencies for 
the spectra and the times for the light curves that 
correspond to the maximum flux calculated numer- 
ically and analytically, in accordance with the above 
estimates. 

The plot of afterglow intensity versus observation 
angle (a ~ 9) (Fig. 6) at different frequencies reveals 
a bright ring attributable to the hotter matter at early 
afterglow stages. The higher the radiation frequency, 
the larger the contrast between the itnage center and 
edge. The same result was obtained by Granot et al. 
(1999). 

The GRB energy must be released in a narrow 
cone, a jet. However, at early stages, the pattern for 
an observer near the cone axis will differ only slightly 
from the pattern produced by a spherical shock. At 
late stages, the jet becomes spherical. 

In Fig. 7, the part of the theoretical light curve 
near the R and K bands is highlighted. We see that 



Fig. 7. Part of the theoretical light curve near the R and 
K bands. 

at longer wavelengths, the light curve passes to a de- 
cline more slowly than it does at shorter wavelengths. 
This chromatic behavior is characteristic of the opti- 
cal afterglows from GRB 990510 (Stanek et al. 1999) 
and GRB 000301c (Jensen et al. 2001 ). 

These objects deserve a more detailed study, but so 
far our model disregards several physical effects (the 
inverse Compton radiation, the Klein— Nishina effect, 
the non-power-law shape of the self-consistent elec- 
tron spectrum, and others). Therefore, it cannot be 
directly used to interpret the spectra of early GRB af- 
terglows. The work to take these effects into account 
is being continued. 

CONCLUSIONS 

The most popular method for analytically estimat- 
ing the spectra and light curves of GRB afterglows 
involves deriving the characteristic frequencies and 
times, determining their behavior, and calculating the 
corresponding flux and constructing the power-law 
segments of the spectra and light curves from the de- 
rived values. These methods are extensively presented 
in the literature (Granot and Sari 2001; Sari et al. 
1998; Wijers and Galama 1999; Waxman 1997). The 
characteristic frequencies in different papers for the 
same cases occasionally differ by a factor of 70, in- 
cluding those in the same papers, suggesting that 
these treatments are insufficient. 

Our calculations unambiguously describe the af- 
terglow spectra and light curves in terms of the model 
under consideration by removing the uncertainty in 
characteristic parameters and without using the next 
approximations. The behavior of our results shows a 
relationship to the observed afterglows, which gives 
confidence that our technique can be used to study 
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the early generation phases of gamma-ray emission 
(during collisions between internal shocks). 
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